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The transport of water in nanoconfined geometries is different from bulk phase and has 
tremendous implications in nanotechnology and biotechnology. Here molecular dynamics is 
used to compute the self-diffusion coefficient D of water within nanopores, around 
nanoparticles, carbon nanotubes and proteins. For almost 60 different cases, D is found to 
scale linearly with the sole parameter 6 as D(6) = D e [1 + (D c /D e - 1)0], with D B and D c the 
bulk and totally confined diffusion of water, respectively. The parameter 9 is primarily 
influenced by geometry and represents the ratio between the confined and total water 
volumes. The D(6) relationship is interpreted within the thermodynamics of supercooled 
water. As an example, such relationship is shown to accurately predict the relaxometric 
response of contrast agents for magnetic resonance imaging. The D(6) relationship can help 
in interpreting the transport of water molecules under nanoconfined conditions and tailoring 
nanostructures with precise modulation of water mobility. 
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Despite its fundamental importance in science and 
technology, the physical and transport properties of water 
are far from being completely understood 1 . The self- 
diffusion of water molecules D in proximity of solid surfaces, 
at the interface between immiscible liquids, and in confined 
geometries, such as nanopores and nanotubes, is a very different 
process as compared to the bulk phase 2-4 . The thermal agitation 
of the water molecules in the bulk liquid is only dictated by the 
local temperature and pressure conditions, and molecular 
diffusion follows the Einstein relation 5 . Differently, under 
confined conditions, the mobility of water molecules is 
perturbed by the presence of additional interaction forces 
arising at the water/solid interfaces, mainly van der Waals and 
Coulomb interactions. These additional forces usually reduce the 
local molecular diffusion 6 ' 7 . Even if considerable work has been 
done in recent years, both experimentally and theoretically, to 
understand and characterize the perturbed behaviour of the water 
molecules in confined geometries, there is still no complete 
comprehension of the process and often the published results are 
contradictory 8 . 

Controlling the mobility of water molecules is of relevance to 
several scientific disciplines and has implications in multiple 
technological applications. For instance, water adsorption/ 
desorption in nanoporous materials, such as zeolites, has potential 
in long-term thermal storage and energy engineering 9 ' ; filters 
with nanopores and nanochannels are increasingly explored for 
their large surface area and higher efficiency 11 ' 12 ; in heat transfer 
problems, nanofluids are under investigation because of their 
peculiar thermal properties 13 ' 14 ; in micro/nanotechnology 
processes, controlling the deposition and surface diffusion of 
water molecules is critical for precise manufacturing 15 ' 16 ; in 
biology, the mechanisms regulating the transport of single water 
molecules through cell membrane channels (aquaporins) and the 
multi- scale water compartmentalization in tissues are still 
elusive 17-20 . Also, proteins tend to modify their structure and 
function according to the surrounding aqueous environment 21 ' 22 . 

Certainly, nanomedicine is one of the fields where several 
exciting discoveries and technological applications can be directly 
related to the anomalous behaviour of water in confined 
geometries. A few examples are the enhancement in longitudinal 
relaxivity associated with the entrapment of Gd 3 + -ion complexes 
in mesoporous structures 23,24 ; the dynamics of water molecules in 
nanotubes and nanochannels for controlled drug delivery 25 ' 26 ; 
and the design of hydrogel-based nano/microparticles 27 ' 28 . 
In particular, the dynamics of water molecules is essential 
in magnetic resonance imaging (MRI), in that contrast 
enhancement is influenced by the local diffusion of water 
molecules 29 ' 30 . It is known that for paramagnetic metal 
complexes, such as Gd 3+ ions, the Solomon-Bloembergen- 
Morgan theory 31 would predict a change in longitudinal 
relaxivity r x of the complex following a variation in the relative 
translational diffusion time (t d ) of the water molecules 
surrounding the complex, and in the residence lifetime (t m ) of 
the water molecules bound to the complex. Similarly, for 
magnetic nanoparticles (NPs), such as the iron oxide NPs, an 
increase in x D (that is, decrease in D) would enhance the 
transversal relaxivity r 2 (ref. 32). Hence, the modulation and 
precise control of the diffusion of the water molecules in the 
vicinity of an MRI contrast agent plays an important role in 
imaging performance. This concept has been already successfully 
proved by experiments 23 , but a clear rationale (and a 
computationally efficient tool) for optimally designing such 
agents is still missing. 

In this work, the self- diffusion coefficient D of water molecules 
is investigated through molecular dynamics (MD) simulations 
under five different isothermal configurations, namely, within 
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silica (Si0 2 ) nanopores, around spherical hydroxylated NPs, 
within Si0 2 nanopores filled by NPs, around single-wall carbon 
nanotubes (CNTs) and proteins. The coefficient D has been 
estimated for almost 60 cases by varying the size of the NPs and 
nanopores, the electrostatic surface charges and level of 
hydrophobicity, as well as the type of protein. The self- diffusion 
coefficient D for all different configurations has been found to 
scale with a single non-dimensional parameter 0, incorporating 
both geometrical and physicochemical information, following the 
relationship D(6) = D B [l + (D C /D B — 1)6]. The D(6) scaling is 
modulated by the coefficients D B and D c , which represent the 
bulk and totally confined diffusion of water, respectively. This 
D(6) law has been applied to estimate the enhancement in MRI 
contrast in magnetic nanoconstructs obtained by geometrically 
confining super-paramagnetic iron oxide NPs (SPIOs) into 
silicon mesoporous matrices. It has been confirmed that the 
transversal relaxivity of SPIOs can be significantly augmented by 
modulating the diffusion of water molecules. This law would help 
in explaining and rationalizing previous experimental evi- 
dences 23 , and represent a ready-to-use tool for the rational 
design of nanoconstructs based on the nanoscale confinement of 
water molecules. 



Results 

Computing the diffusion of nanoconfined water molecules. 

MD simulations were used to compute the self- diffusion 
coefficient D of water molecules confined under different con- 
figurations. These are shown in Fig. 1 and include the case of 
water molecules (blue dots) moving (a) around spherical hydro- 
xylated nanoparticles (NPs) (grey dots); (b) within a hydrated 
nanopore (grey dots); (c) around hydroxylated NPs (red dots) 
adsorbed on the surface of a hydrated nanopore (grey dots); (d) 
around and within single- walled carbon nanotubes (CNTs); (e,f) 
around proteins. The NPs are made out of magnetite (Fe 3 0 4 ) 
crystals (red and cyan dots), with OH - functional groups on 
their surface, or Si0 2 crystals (grey dots), with silanol SiOH 
functional groups on the surface. The nanopores are made out of 
Si0 2 only. 

To investigate the influence of geometry and material proper- 
ties, the self- diffusion coefficient D of the water molecules was 
computed for 58 different cases. In particular, these cases were 
different in terms of NP diameter, being 1.3, 2.0 or 5.2 nm; 
nanopore diameter, with values 2.0, 4.1, 8.1 and ll.Onm; number 
of NPs adsorbed on the nanopore wall, varying from 0 to 66 NPs 
per pore; CNT armchair chirality, namely, (5,5), (10,10), (20,20) 
and (30,30); and type of proteins, including molecules with a 
spherical (for example, ubiquitin) and elongated (for example, 
Ca 2 + -ATPase) shapes; of small (for example, 562 atoms of 
Bl-immunoglobulin-binding domain) and large (for example, 
9667 atoms of Ca 2 + -ATPase) sizes; and exhibiting a catalytic 
(for example, glucokinase), hormonal (for example, leptin) and 
transport (for example, myoglobin) function. As per the material 
properties, two types of interactions were considered in the MD 
simulations: bonded interactions and non-bonded interactions 
between the water molecules and the solid surfaces, described via 
van der Waals and Coulomb potentials. In the Si0 2 structures 
and spherical NPs, the bonded interactions are modelled by 
means of harmonic terms. For the NPs, the strength & of the 
Lennard- Jones potential was varied from 2.49 to 24.94 kj mol _ 1 
and the partial electrostatic charges of atoms were set to either the 
nominal value or zero for NPs and nanopores. For simulations 
with CNTs, the non-bonded interactions between CNTs and 
water molecules were modelled by the Lennard-Jones potential, 
with neutral carbon atoms and o cc — 0.36 nm, e cc = 0.29 
kjmol -1 (ref. 33). Finally for the proteins, all bonded and 
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non-bonded interactions were modelled using the GROMOS96 
43a2 force field 34 , which has been widely used for studying 
similar applications 35 ' 36 . Finally, various hydration levels were 
considered providing an overall water density ranging from 715 
to 941 kg m -3 . Note that for these density values, there is no 
heterogenous wetting and consequently no anomalous behaviour 
related to low water-filling regimes 3 ' 7 . 

All computed values of the self- diffusion coefficient D are 
reported in the Supplementary Tables 1-11 and Supplementary 
Note 1. In general, it is observed that the coefficient D decreases 
compared with the bulk value, 2.60 x 10 _9 m 2 s~ 1 at 300 K 
(ref. 37), as the ratio between the total area of the solid-liquid 
interface and the total volume occupied by the water (V w ) 
increases. More specifically, for a Si0 2 nanopore, D is reduced 
from 2.50 ± 0.09 x 10- 9 m 2 s _1 to 0.82±0.22 x 10- 9 m 2 s _1 as 
the pore diameter decreases from 11.0 to 2.0 nm. Moreover, D is 
inversely proportional to the NP concentration and diameter. For 
a 5.2-nm Si0 2 NP, D decreases from the almost unconfined value 
to 2.12 ± 0.04 x 10 _9 m 2 s _1 following a 36% increase in NP 
concentration. Indeed, the increase in NP concentration is 
associated with a decrease in separation distance between 
adjacent NPs and, consequently, a decrease in the volume 
available to water molecules. Consistently, D reduces as the 
concentration of Fe 3 0 4 NP increases within a nanopore: D is 
2.2010.10 x 10 _9 m 2 s _1 in an 8.1-nm Si0 2 nanopore; how- 
ever, if 2 Fe 3 0 4 NPs, of 2.0 nm in diameter, are adsorbed on the 
nanopore surface, D decreases to 2.07 ± 0.14 x 10 _ 9 m 2 s _ 1 and 
it drops to 0.44 ± 0.05 x 10~ 9 m 2 s _1 (-80% decrease), if 16 
Fe 3 0 4 NPs, of 2.0 nm in diameter, are added in the nanopore. In 
addition, for 16 Fe 3 0 4 NPs adsorbed on the wall of an 8.1-nm 
nanopore, D is 1.46 ± 0.09 x 10 _9 m 2 s _1 for 1.3 nm NPs and 
becomes 0.44 ± 0.05 x 10 ~ 9 m 2 s " 1 ( ~ 70% decrease) for 2.0 nm 
NPs. A similar trend is observed for the CNTs and proteins by 
reducing the size of the water box. More specifically, around the 



(5,5) chirality CNT with a length of 5nm, the water diffusivity 
D decreases from the bulk value (box of 316 nm 3 ) to 
1.22 ± 0.12 x 10~ 9 m 2 s _1 (box of 21 nm ). Similarly, around 
the Bl-immunoglobulin-binding domain, the diffusivity D 
decreases from 2.41 ± 0.04 x 10 _9 m 2 s _1 (box of 348 nm 3 ) to 
0.87 ± 0.10 x 10 ~ 9 m 2 s _1 (box of 23 nm 3 ). As expected, the 
above data qualitatively demonstrate that the self- diffusion 
coefficient D of water is strongly correlated to the ratio between 
the interface surface and the total water volume: the larger this 
ratio the smaller the water mobility. However, the possible 
contribution of other parameters should also be assessed. 

To this end, sensitivity analyses were performed to elucidate 
the effect of the Lennard- Jones potential strength s and Coulomb 
interactions on D. Larger values of the parameter s are associated 
with lower mobility of the water molecules. As an example, let us 
consider the case of eight Fe 3 0 4 NPs (2.0 nm diameter) adsorbed 
on the walls of an 8.1-nm diameter Si0 2 nanopore. A one order 
of magnitude decrease of a of Fe 3 0 4 atoms only carries a 10% 
increase of D: for s = 24.94 kjmol - 1 D is equal to 1.33 ± 0.13 x 
10 ~ 9 m 2 s ~ l , whereas D increases to 1.47 ± 0.1 1 x 10 ~ 9 m 2 s ~ 1 
for £ = 2.49 kj mol - 1 . Moreover, D increases as the surface 
electrostatic charges decrease. For neutral Fe 3 0 4 NPs, D grows to 
1.64 ± 0.04 x 10 _9 m 2 s _1 and, if both the NPs and nanopore 
wall are electrically neutral, D takes the value of 1.69 ±0.20 x 
10 _9 m 2 s _1 ( — 30% increase as compared with the example 
above). Although the water molecule confinement is affected by 
the strength of the interaction potentials (van der Waals and 
Coulomb), geometrical parameters show a greater influence on the 
coefficient D. The reason is that all considered surfaces have effective 
wall potentials that are strong enough to induce a significant 
reduction of the water mobility in a region close to the wall. On the 
other hand, as clarified below, the volume of the low mobility region 
only slightly depends on the wall potential strength, namely, the 
minimum of the potential well generated by the wall. 



a b c 




Figure 1 | Selected configurations, (a) Si0 2 particle in water, diameter </> = 5.2nm (blue dots: water molecules; grey dots: Si0 2 atoms); (b) Si0 2 
nanopore filled by water, diameter <E = 8.1nm; (c) sixteen Fe 3 0 4 NPs within a Si0 2 nanopore filled by water, cf) = 2.0nm and <I> = 8.1nm (red and cyan 
dots: Fe 3 0 4 atoms); (d) single-walled CNT with chirality (5,5); (e) green fluorescence protein; (f) leptin protein (the standard ribbon visualization of 
secondary structures has been used for proteins). In d-f water molecules have been removed for clarity. Almost 60 different cases have been analysed 
by varying the size and surface properties of the NPs, nanopores and nanotubes as well as the type of protein. 
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Finally, the self- diffusion coefficient D did not change 
significantly with the level of hydration in the considered range, 
in accordance with previous studies 7 . Considering again the 
representative case of eight Fe 3 0 4 NPs, of 2.0 nm in diameter, 
adsorbed on the wall of an 8.1-nm Si0 2 nanopore, D ranges 
from 1.30 ± 0.11 x 10~ 9 to 1.40 ± 0.07 x 10 " 9 m 2 s " 1 (<10% 
variation) as the water density increases from 700 to 930 kg m -3 . 

Characteristic length of confinement. In the bulk fluid, the water 
molecules fluctuate with a kinetic energy proportional to k B T y 
where k B is the Boltzmann constant (1.38 x 10 _23 JK _1 ) and 
T is the temperature. As opposed to molecules in the bulk, those 
in a close proximity of solid surfaces are subjected to additional 
van der Waals (C/ V dw) an d Coulomb (U c ) interactions interfering 
with their state of agitation. This induces a layering of water 
molecules with reduced mobility near the solid surface (see 
Supplementary Figs 1-3 and Supplementary Discussion), as 
already pointed out in other works ' 39 . A characteristic length 
S can be introduced to quantify the thickness of such confined 
water layer. 

Referring to the popular notion of solvent accessible surfaces 
(SAS) 40-42 , the quantities S tot and Si oc can be introduced as the 
total and specific (per atom) SAS areas, respectively. For an 
arbitrary atom i of the solid structure, a number N n of nearest 
neighbours (including the atom / itself) can be identified within a 



fixed cutoff radius (Fig. 2a,b). The corresponding effective 
potential energy (7 eff on the water molecules, due to both van 
der Waals (l/ v dw) and Coulomb (U c ) interactions, can be 
computed as: 

Utf(n) = U y6w (n) + (U c )(n), (1) 

along the n direction, orthogonally to the SAS and passing 
through the centre of the atom i (Fig. 2b). For the 12-6 
Lennard- Jones potential, it follows that U v d w (n)= Ylk=i ^ s k 

(nf) ~ (?) ' wit ^ Sh ° k anc * Tk d enotin & tne depth of the 

potential well, the distance at with such potential becomes zero 
and the Euclidean distance between the generic line point with 
coordinate n and the centre of kth. nearest neighbour, respectively. 
For the Coulomb interactions, the average potential energy at a 
fixed temperature T between the N n atoms and the water 
dipoles is (U c )(n) = — Efi w r(Efj, w /k B T) y where £, /i w , k B and r 
denote the electrical field strength, water dipole moment 
(7.50 x 10 ~ 30 Cm for the SPC/E model), the Boltzmann 
constant and the Langevin function r(x) — coth(x) — 1/x. The 
strength of the electrical field E can be readily computed following 
the law of electrostatics (see Methods). Knowing the effective 
potential U e f£n) for the atom /, a corresponding characteristic 
length Si can be estimated within which the water molecules have 
reduced mobility. This length Si is given by S t — n i>2 — «f,i where 




Figure 2 | The characteristic length of confinement 3. (a) A local characteristic length 3 } can be defined at any atom /' (with non-vanishing SAS). 

(b) The contribution to the total potential energy of neighbours along a direction orthogonal to the SAS (line / with local coordinate n) is computed. 

(c) The thermal energy level provides a criterion to (locally) define the characteristic length 3„ which typically varies along the whole SAS. 

(d) A characteristic length of confinement 3 for the whole structure can be defined by a weighted average of all 5,. (e) The length of confinement is 
reported for a Si0 2 nanopore (^^0.33nm) and for a Fe 3 0 4 NP ((5 2 ~0.50 nm). For Si0 2 , only the first water layer (with distance r w1 from the pore axis) 
is located within the volume of influence. For Fe 3 0 4 , both the first and the second water layers (with distances r w1 and r w2 from the particle centre, 
respectively) are found within the volume of influence (more cases are reported in the Supplementary Figs 1 and 2). 
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n^ 2 and n iy i are the two zeros of the equation U e ^n) + k B T/4 = 0 
(Fig. 2c). Therefore, based on the definition of S h all the water 
molecules located within such a distance are significantly affected 
by the van der Walls and Coulomb interactions, whereas all the 
water molecules beyond the characteristic length Si can escape the 
potential well generated by the solid wall. Here k B T/4 is half the 
kinetic energy per independent degree of freedom according to 
the equipartition of energy (see Methods). By proper averaging 
over the surface, the mean characteristic length S of the overall 
solid surface (Fig. 2d) can be derived as 



N 

<5;Sioc,z 
.i=i 

Stot 



(2) 



with Si OCj ; and N being the specific (per atom) SAS for the atom 
i and the total number of atoms, respectively. Note that the 
above formulation is general and applies to hydrophilic and 
hydrophobic surfaces, regardless of their electrostatic surface 
charge. Also, the characteristic length S can be conveniently 
computed based on the geometry of the problem, Lennard- Jones 
force field parameters and the partial electrostatic charges using 
the script provided in the Supplementary Software 1. 

In Fig. 2e, the water density profile within (a) a Si0 2 nanopore 
and (b) around a single Fe 3 0 4 NP is shown, where peaks denote 
the typical water layers nearby a solid wall at the nanoscale 2 ' 39 . 
It is noteworthy that despite the great difference in the potential 
strength between Si0 2 and Fe 3 0 4 (min((7 e ff j2 )/min((7 e ff )1 ) «2.7), 
the corresponding difference in terms of characteristic lengths is 
much more moderate (<5 2 A5i^ 1-5). From the density profiles, 
Fe 3 0 4 clearly induces a stronger perturbation in the nearby water 
molecules' distribution. However, except for a first thin water 
layer strongly adsorbed to the NP surface (and accounted by 
S 2 >Si) y the amplitude of these perturbations rapidly decays 
further away and becomes comparable in the remaining confined 
volume for both cases. Similarly, although there is a significant 
difference in the potential minimum between green fluorescence 
protein and (5,5) chirality CNT (min(L7 e ff >2 )/min(L7 e ff >1 ) « 1.5), the 
difference between the two characteristic lengths is negligible. The 
above observations suggest that geometrical parameters could 
more significantly affect D as compared to energetic parameters. 

Scaling law. Since water mobility is impaired mostly in a thin 
layer next to the liquid- solid interface with thickness d, it is 
reasonable to assume that the observed variation in the self- 
diffusion coefficient D is mainly associated with the altered 
mobility of the water molecules within such a layer and the 
corresponding volume (volume of influence — Fig. 2d). In addi- 
tion, it has been already observed that the self- diffusion coeffi- 
cient D reduces as the ratio between the total interfacial area and 
the total volume occupied by water increases. On the basis of such 
an evidence, a scaling parameter 6 can be introduced as the ratio 
between the total water volume of influence (V in ) and the total 
volume accessible to the water molecules (V w ); thus, 



Vin p 



(3) 



This parameter 6 varies from 0 (bulk water case) to 1 (totally 
confined water). The volume of influence (y in ) is the volume 
of water that feels the van der Walls and Coulomb interactions 
and is therefore influenced by the presence of solid walls. This 
volume is readily given 

by V in Kj2 p s{P>5(P\ where S$ and 5® 
represent the total SAS and characteristic length of the pth 
particle, respectively. As detailed in the Methods, possible overlap 
of the volumes of influence due to several particles/pores (for 



example, several NPs loaded within a Si0 2 nanopore) can be 
easily taken into account by the continuum percolation theory 
(CPT), and a more accurate estimate of Vi n may include the effect 
of particle curvature. 

Assuming 6 as the sole, independent variable for D, all 
computed values relax within a narrow band around a linear 
curve (Fig. 3) that can be readily described by the relationship 



D(ey 



(4) 



where D B is the self- diffusion coefficient of bulk water, while D c 
the self- diffusion coefficient of totally confined water. Remark- 
ably, Fig. 3 presents data from 58 different cases analysed in this 
work as well as data available in the published literature. Here 
despite the variety of the considered configurations, particles and 
sources of the results, a simple law is found to be sufficiently 
accurate to describe the phenomenon under study, thus 
confirming that 6 is indeed an important controlling parameter 
under very diverse conditions and geometrical configurations. To 
have a more explicit formulation of equation (4), details about the 
evaluation of D c are provided in the section below. 

Thermodynamic insights. Nanoconfined water shares some 
features with supercooled ordinary water in that it may not 
crystallize on cooling below the melting temperature of 
T M ^273.15K (refs 1,43). Within the thin S layer of water 
molecules next to a solid surface, the thermodynamic state 
depends on the characteristic confinement length scale 44 . In 
particular, the specific heat capacity c p of nanoconfined water has 
been experimentally measured in narrow Si0 2 nanopores and its 
variation with the temperature T is plotted in Fig. 4a for a pore 
diameter of 1.7 nm (ref. 45). Using these experimental data, the 
energy variation associated with the transition from bulk 
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Figure 3 | Scaling behaviour of the water diffusion coefficient D. The 

self-diffusion coefficient of water D has been calculated for over 60 
different cases including spherical NPs in water, water within nanopores 
with and without spherical NPs, proteins and CNTs in water. Fifty-eight 
cases are the results of MD simulations performed in this work, under 
isothermal conditions. Results from the literature are also provided, for 
which the scaling variable 6 was computed as from the Supplementary 
Discussion. The solid and the dashed lines represent equation (4) for 
D c = 0 and D c = 0.39 x 10 ~ 9 m 2 s ~ \ respectively. Equation (4) accurately 
recovers simulation and literature results with high coefficient of 
determination (ft 2 >0.90). The uncertainties on the value of D (vertical 
bars) refer to the fitting of the mean square displacement, whereas the 
uncertainties on the value of 6 (horizontal bars) refer to the estimate of 
the total volume accessible to water molecules V w (see Methods). 
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Figure 4 | Derivation of the function D C /D B = g( - s). (a) The specific heat capacity c p (T) of water confined in the narrowest Si0 2 pore reported in Nagoe 
et a/. 45 (b) By integrating the latter function c p (T), the energy of nanoconfinement is eva luated as Ah=f(T)=- f T ° c p 6T. (c) The inverse function 
T = f _1 (A/i) provides the temperature corresponding to a given energy Ah; hence, it can be used in combination with the experimental results from 
Chen et a/. 1 where the diffusion coefficient (in narrow Si0 2 pores) is given as a function of T. (d) The ratio between the diffusion D c of fully confined 
water and the bulk diffusion D b of water, D C /D B = g(Ah), is reported and the corresponding values for all the considered configurations are shown, 
where g(A/i)<0.15. 



to confined water can be readily computed as 
Ah=- f T ° CpdT=f(T), with /(T o ) = 0 at the bulk water 
temperature T 0 = 300K. From this, the inverse function 
T—f~ l {Sh) can be derived as plotted in Fig. 4b. The energy 
variation associated with water nanoconfinement can also be 
written as Mi— he — h B <0, where h B and h c are the energy of 
bulk and confined water, respectively. Since the confined region is 
typically limited to 1-2 layers of water molecules (for all the 
considered structures, S < 0.6 nm), it can be assumed that h « — e, 
which is the minimum of the effective potential (l/ eff ) generated 
by the solid surface. The function/ - 1 can be used to compute the 
supercooled temperature of nanoconfined water corresponding to 
the energy variation — e, that is, to say T—f~ 1 ( — s). In the work 
by Chen et al. 1 , the diffusion coefficient of water molecules under 
strong nanoconfinement (that is, within Si0 2 pores with diameter 
of 1.4-1.8 nm) is reported as a function of the temperature. From 
the latter experimental data, the diffusion coefficient of totally 
confined water D c can be readily expressed in the form 
D C — D B g(T) (Fig. 4c). Therefore, by combining the two sets of 
experimental data and knowing the value of e, the diffusion 
coefficient of totally confined water D c can be computed as 
D c /D B =g(T)=g(f- 1 (-8))=g(-s). The ratio D C /D B is 
plotted in Fig. 4d for the different 58 cases analysed here. For 
the iron oxide NPs and CNTs, which are characterized by a fairly 
strong effective potential (large e), D c ^0. As the strength of the 
effective potential reduces, the value for the diffusion coefficient 
of totally confined water increases. However, even for the Si0 2 
NPs and nanopores, and for the proteins that are characterized by 
lower e, in our computations D c is at most 15% of D B . We stress 
that the function g(Ah) is a property of water, and here we have 
chosen to estimate it on the basis of measurements in Si0 2 
nanopores only, because such experiments are among the very 
few that are well documented in literature. 

In general, the volume of water can be partitioned into bulk (B) 
and confined (C). Invoking the mixing rule, the average 
diffusivity D of the system can be presented as 



where N c and N B are the number of water molecules in the 
confined and bulk regions, respectively. Considering that AT C = 
p c y c /m w , with p c being the water (mean) density in the 
adsorbed region and ra w being the mass of one water molecule, it 
follows: 



D(6,s)=D B \ 1 + 



Dc 



1 



PcO 

Pc0 + Pb(1-0), 



(6) 



D= 1 



N c 



N B +N C 



D B + 



N c 



N B +N C 



-Dc 



(5) 



This implies that the average diffusivity D depends in general 
on a geometrical parameter (9) and an energetic parameter (s). 
However, the following approximation (pc&)Kpc& + Pb 
(\ — 6))&6 can be safely made for the full range of 6 (see 
Supplementary Fig. 4 and Supplementary Discussion). Therefore, 
equation (6) degenerates into equation (4) demonstrating that the 
diffusion of nanoconfined water can be interpreted invoking the 
thermodynamics of supercooled water. Moreover, based on the 
above discussion on D c (see Fig. 4d), for a large variety of 
nanoconfined systems, the simplifying assumption D C /D B « 0 can 
be safely made. In such cases, the much simpler law 
D(6,&)&D(6) = D B (1 — 0) is readily derived showing the direct 
linear dependence of the water diffusion coefficient on the sole 
scaling parameter 6. Note that no empirical factor is needed to 
derive D(9) y with the latter law matching the values of the 58 MD 
simulations as well as 13 further configurations from the 
literature with a quite good coefficient of determination 
(.R 2 = 0.93, solid line in Fig. 3). The dashed line in Fig. 3 
represents, instead, equation (4) with D C ID B — 0.1 5, correspond- 
ing to the largest value of D c only observed in a few simulated 
cases (Fig. 4d). 

Discussion 

It has been shown that the relaxometric properties of contrast 
agents for MR imaging can be enhanced by proper modulation 
of the water molecule dynamics 46 . Recently, Decuzzi and 
collaborators 23,24 have shown that the geometrical confinement 
of Gd 3 + -based MRI contrast agents and SPIOs within 
mesoporous structures can increase by several folds the 
longitudinal and transversal relaxivities, respectively, of the 
original agents. Here it is shown that the scaling law D(9) can 
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Figure 5 | Relaxivity of SPIO-loaded mesoporous silicon nanoconstructs. (a) Scanning electron microscopy image of discoidal mesoporous silicon 
nanoconstructs (with cross section area 1000 x 400 nm 2 ) loaded with 5 nm SPIOs. Scale bar, 1 |im. (b) Schematic representation of the pores within the 
silicon nanoconstructs loaded with the SPIOs (red circles), (c) TEM and EDX images showing the distribution of SPIO clusters within the silicon 
nanoconstructs. The black dots in the TEM image and red dots in the EDX image identify the SPIOs; the yellow lines in the EDX image represent the Si walls 
of the mesoporous silicon nanoconstructs. Scale bar, 100 nm. (d) Diagram showing the local enhancement in transversal magnetic resonance relaxivity as 
calculated via the D(6) relationship. 



be effectively used to predict the enhancement in transversal 
relaxivity exhibited by SPIOs confined in mesoporous silicon 
structures. Referring to Fig. 5a, discoidal mesoporous silicon 
particles (SiMPs) of 1,000 nm in diameter and with an average 
pore size of 40 nm are synthesized by a combination of optical 
lithographic techniques and electrochemical etching 26 . These 
particles are exposed to a concentrated solution of SPIOs, which 
are loaded within the pores by capillary suction and form clusters 
of controlled size. As schematically shown in Fig. 5b and 
documented by the transmission electron microscopy (TEM) and 
energy-dispersive X-ray spectroscopy (EDX) images of Fig. 5c, the 
SPIOs are distributed more or less uniformly within the pores and 
across the porous matrix of the silicon particles. The black and 
red spots, respectively, in the TEM and EDX images identify the 
SPIO clusters within the mesoporous silicon matrix. In these 
experiments, commercially available SPIOs (Sigma-Aldrich) are 
considered, coated with a thin layer of polyethylene glycol to 
facilitate their dispersion in water, and presenting an average core 
diameter of 5.13 ± 1.07 nm, as derived by TEM analysis. In bulk 
water the transversal relaxivity of the unconfined SPIOs has been 
measured to be (r 2 = 107 ± 24 mM ~ 1 s _ whereas on geometri- 
cal confinement within the SiMPs, the relaxivity rises up to 
r 2 = 270±73 mM _1 s _1 . 

From the theory on the MR relaxation of iron oxide NPs, the 
transversal relaxivity r 2 of the SPIOs dispersed in an aqueous 
solution can be estimated as 29 ' 32 ' 47 



r 2 = 



Mi 



i 

To* 
Fe 



■i. 1 V 

with — oc — , 
T 2 D 



(?) 



where T 2 and T 0jW are the transversal relaxation times for the 
solution with contrast agents and bulk water ( = 2,800 ms), 
respectively; M Fe is the iron concentration (mM) in the solution; 
v is the volume fraction of the SPIOs; and D is the diffusion of the 
water molecules surrounding the SPIOs. Details on equation (7) 
are provided in the Supplementary Discussion. 

For a fixed iron concentration M Fe in the solution, the 
enhancement in transversal relaxivity is directly related to the 
ratio En = (v/D)/(u B /D B ), where the quantity v/D has to be 
estimated for confined SPIOs and v B /D B for the free SPIOs 
dispersed in bulk water. Since the distribution of the SPIOs within 
the SiMP is not uniform, the diffusion coefficient D and the local 
SPIO volume fraction v are expected to vary within the porous 



matrix. The local enhancement in relaxivity En is plotted in 
Fig. 5d. Areas are shown with mild and high enhancement, up to 
20-fold, corresponding to different levels of SPIO loading and 
clustering, and therefore different levels of water confinement. By 
integrating over the whole silicon particle, an average relaxivity 
enhancement En can be computed as being ~2.7, which is in 
excellent agreement with the experimental value of ~2.52 (see 
Methods for details). Note that the proposed law can be used 
in several other fields to analyse different cases of water 
confinement. 

In summary, it has been shown that the self-diffusion 
coefficient of nanoconfined water can be described by a unique 
dimensionless parameter 9 y representing the ratio between the 
confined and total water volumes. The coefficient D scales linearly 
with 6 and can be readily estimated, knowing the bulk D B and 
totally confined D c diffusion of water. This has been validated on 
the basis of almost 60 different cases and 5 different geometrical 
configurations, including the analysis of the water molecule 
dynamics within nanopores and CNTs, around NPs and proteins. 
The coefficient of diffusion for confined D c water is quantified on 
the basis of the thermodynamics of supercooled water. As an 
example, the scaling relation has been shown to accurately predict 
the enhancement in magnetic resonance relaxivity of iron oxide 
NPs confined into mesoporous silicon structures. 

The proposed approach may be used to interpret experimental 
data collected in different scientific disciplines on the dynamics of 
water molecules under confined conditions and to rationally 
design nanostructures for modulating the diffusion of water. This 
is of relevance in nanomedicine, nanotechnology as well as in 
more traditional engineering fields such as heat transport, fluid 
dynamics and energy storage. 

Methods 

MD simulations. Atomic coordinates of Fe 3 0 4 NPs are generated from Fe 3 0 4 
crystals 48 , whereas Si0 2 crystals are considered in the case of nanopores and Si0 2 
NPs (Supplementary Figs 5-8) 49 . Crystal structures of proteins are taken from the 
Protein Data Bank (http://www.rcsb.org; Supplementary Figs 9 and 10). CNTs are 
generated by means of the Visual Molecular Dynamics software (Supplementary 
Fig. II) 50 . The SPC/E model 51 is used for water molecules, which is known to 
accurately predict some of the properties of water relevant for this study at room 
temperature 37 . However, it is also worth noticing that the SPC/E model does not 
accurately predict some other properties of water. For instance, shear viscosity or 
thermal conductivity were found to be off by more than 50% at room 
temperature 52 . Bonded interactions of Si0 2 and Fe 3 0 4 are treated by means of 
harmonic stretching and angle potentials 53 . Van der Waals interactions are 
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modelled by a 12-6 Lennard- Jones potential; partial charge interactions between 
solid surfaces and water are modelled by a Coulomb potential (Supplementary 
Tables 12- 14) 51 ' 54-56 . Non-zero partial charges are only assigned to atoms on the 
surface of nanopore or NP, which belong to silanol and FeOH groups, whereas all 
other atoms (bulk of Si0 2 and Fe 3 0 4 ) are considered as neutral. 

Simulations are carried out with a leap-frog algorithm (time step: Af = 0.5fs), 
and periodic boundary conditions are applied along the three Cartesian 
coordinates. After energy minimization of NP, nanopore or nanotube setups, the 
two subsystems (solid crystals and water) are initialized at 300 K (Maxwellian 
distribution of velocities) and fully coupled to a Nose-Hoover thermostat 57 ' 58 
(at 300 K and time constant x = 0.2 ps) for 50 ps, until the energies of the system 
relax to a steady state. During the latter preliminary calculation, one thermostat for 
each subsystem is adopted. Afterwards, Nose-Hoover thermostats (at 300 K) are 
maintained attached to solid crystals only, whereas the simulation is continued up 
to 2 ns. In the case of proteins, the MD protocol is slightly changed to improve 
convergence. Energy minimization for proteins is performed before and after 
solvation and ions are added when needed for achieving the neutrality of the 
system, which is then equilibrated in two steps: 100 ps in canonical ensemble (fixed 
number N of particles, volume V and temperature T - NVT) at 300 K (initialization 
with Maxwellian distribution of velocities, Berendsen thermostats 59 separately 
attached to proteins and water, x = 0.1 ps); 100 ps in NPT ensemble (fixed number 
N of particles, pressure P and temperature T) at 300 K and 1 bar (Berendsen 
thermostats separately attached to proteins and water, t = 0.1 ps; Parrinello- 
Rahman pressostat 60 applied to the whole system, T = 2ps). During the 
equilibration, all bonds in the proteins are kept rigid using the LINCS (Linear 
Constraint Solver) algorithm 61 . Finally, a Nose-Hoover thermostat (300 K, 
x = 0.2 ps) is attached to protein's atoms and the simulation is continued for 1 ns. 

In all simulated cases, steady state is reached when D, which is evaluated every 
200 ps, tends to an asymptotic value (Supplementary Figs 12-28 and 
Supplementary Note 2). This is generally achieved after about 600 ps, for all 
configurations. Note that the root mean square deviation of the proteins in the 
water box with respect to the crystallographic structures is on average found to be 
below 0.3 nm (Supplementary Fig. 29 and Supplementary Note 2). The self- 
diffusion coefficient D of the water molecules is determined following the classical 
relationship of Einstein and computing the mean square displacement as 5 ' 62 
MSD = lim ( ||rj(f) — r ; (0) || 2 ) iewater =6i)f, where the position vector refers to the 

centre of mass of the water molecule i at the generic time t and 0 refers to the initial 
configuration of the system. Alternative approaches for computing the water 
diffusivity could be considered as well (for example, those based on the 
first-passage concept) 63 ' 64 . Further details on the implementation of the MD 
simulations and calculation of the self- diffusion coefficient D are provided in the 
Supplementary Methods. 

The MD simulations are performed with the software package GROMACS 65 ' 66 . 
Rendering is performed with UCSF Chimera 67 . 

Derivation of the characteristic length d. When denning the van der Waals 
potential U vdw (n) in equation (1), the parameters s k and o k already incorporate a 
combination rule for the Lennard-Jones parameters between the atom i and 
oxygen atoms of water (for example, the Lorentz-Berthelot rule). Moreover, in 
the average potential (U c )(n)= — £ i w w r(£/i w //c B T), the effective strength of the 
electrical field E= ||£(n)|| may be expressed by the following explicit form: 



E(n) 



1 



47l£ 0 £r 



r k (n] 



r k (n)=^ [x(n) -x k } 2 + \y(n) -y k ] 2 + [z(n) - z k ] 2 



(8) 



(9) 



with q k being the electric charge of the kth neighbour, while (x,y,z) and (x k ,y k ,z k ) 
represent the Cartesian coordinates of the generic point on / (corresponding to the 
local coordinate n) and the Cartesian coordinates of neighbour k, respectively. In 
our approach, the relative permittivity e r is an input parameter to be provided to 
the Matlab(r) routine for the computation of 3 (see below, Supplementary Fig. 30, 
Supplementary Methods and Supplementary Software 1). In particular, in this 
work, e r was included as a function of the distance from the particle £ r (n) following 
the suggestion in ref. 68. 

The expression of (U c )(n) is a classical result of electrostatics and can be 
justified as follows. Let E be the electric field acting on a dipole ji. The 
instantaneous energy can be expressed as U c ((p)= — \i ■ E= — piE cos cp, with E and 
fi being the field and dipole strength, respectively, while q> is the angle between the 
direction of the dipole and the field. The minimum and maximum values of energy 
are attained for (p = 0 and n, respectively. However, in the presence of thermal 
agitation, the direction of \i (hence (p) is continuously changing in time. For 
classical systems, a number of independent dipoles are distributed according to 
their energy level U c :U Cymin < U c <U Cjmax . In particular, at thermodynamic 
equilibrium, the Boltzmann distribution predicts that the number N(U C ) of dipoles 
with energy U c is: N(l/ c )=Cexp( — ^), with C being a constant to be determined. 
Note that, in three dimensions, all dipoles lying on a cone with angle 2<p around the 
electrical field direction have the same energy. Moreover, at every cp, the total 
number of such dipoles is N(U c ((p))d£l, with Q being the incremental solid angle. 



The above dipoles at cp give the component sum d/z in the field direction: 

dpi = NdQfi cos (p. Hence, according to the Boltzmann distribution, the average 

dipole moment takes the more explicit form: 



Jexp(^^)sin(pcos(pd(p 



Jexp(^)sin<pd<p 



(10) 



On substituting fi = piE/k B T and q= cos (p, the above integral becomes 
(n)=liT(P). 

On construction of the potential in equation (1) at an arbitrary atom i, a local 
characteristic length 3 t can be denned as 3 { = n i>2 — n i}1 , where n i>2 and n iA are the 
two zeros of the equation (see Fig. 2c): t/ eff (n) + txk B T= 0, where we expect a« 1/4. 
In fact, provided that k B T/2 is the kinetic energy attributed the each degree of 
freedom of the water molecules, for planar surfaces a = 1/4 because molecule are 
allowed to escape the potential well only along half of the direction orthogonal to 
the surface. The previous equation implies that water molecules located beyond the 
characteristic length 3 { are in the position of escaping the potential well generated 
by the N n atoms in the solid wall owing to the kinetic energy k B T. Obviously, when 
the horizontal line U-\-uk B T = 0 does not intersect the function in equation (1), 
<5,- = 0. 

In general, the quantity 3 t varies at each atom i (see Fig. 2a). Moreover, 
(meaningless) non-zero values for 3i can be found for bulk atoms. Thus, for a given 
solid structure, it is convenient to define a mean characteristic length as in 
equation (2) (see Fig. 2d). Note that both S tot and S\ OCji are readily computed by 
GROMACS, once the geometry of the system is known (for example, in the form of 
a pdb file). It is also worth emphasizing that 3 is a characteristic length of the whole 
system of interest and can be straightforwardly computed based on the geometry, 
Lennard-Jones force field parameters and partial charges. 

The script for computing 3 is based on Matlab(r) and it is provided as 
Supplementary Software 1 of this work. 

Derivation of the scaling parameter 0. Finding a proper scaling parameter 6 
(or equivalently 3) is not trivial. A few unsuccessful attempts to find a general 
scaling parameter for water self- diffusion coefficient are reported in Supplementary 
Discussion and Supplementary Figs 31 and 32. 

In equation (3) the particle curvature is neglected; thus, this is accurate in the 
limit 3^ « with being a representative radius of curvature of the pth 
particle. If the latter approximation is not properly fulfilled, a more accurate 
estimate of the numerator at the right-hand side of equation (3) can be easily 
adopted (see also a few examples in the Supplementary Discussion). 

Clearly, in cases of strong confinement (and high values of 6) with the presence 
of several particles (for example, the reported studies where a number of spherical 
NPs are loaded within cylindrical nanopores), a partial overlap of the volumes of 
influence becomes more and more probable. However, we notice that volumes of 
influence due to different particles are not additive. As a result, the quantity 6 
computed by equation (3) is only apparent, with the effective fraction of the 
volume of influence being smaller than 9. The above issue is encountered in the 
framework of CPT. To this respect, a classical result of CPT, under the assumption 
of randomly placed volumes suggests that, to properly recover the effective fraction, 
the apparent volume fraction should be corrected as: 1 — exp ( — #) 69 ' 70 . Hence, we 
suggest that the above correction applies in the presence of overlap of the volumes 
of influence (for example, several particles within the computational box). 

The volume V w . For fully hydrated simple geometries, V w can be computed by 
considering the nominal size of a particle/pore (see the discussion below about 
relaxivity enhancement). However, in general, the volume V w in equation (3) can 
be denned as: 



V w =N sol /p n , 



(11) 



where iV sol and p n are the number of water molecules within the computational 
periodic box and the average water density, respectively. As a consequence, the 
volume occupied by a solvated particle p is V p = V hox — V w , where V box is the 
volume of the computational box. For complex configurations (for example, several 
NPs surrounded by water or within nanopores), the volume V w can be estimated as 
^w=^out — J2pL\ Vp'K where and N p are the volume of the pth particle and 
the total number of particles, respectively, whereas V out is the volume of the 
surrounding space (that is, y out = V box for particles not loaded in nanopores, while 
^out = Vp 0re f° r cases where particles are loaded within a pore whose volume, 
computed by equation (11), is Vp 0r e)- 

Clearly, the use of equation (11) relies on the computation of the number 
density p n . Using available packages in standard MD software, an estimate of the 
average value for p n can be easily computed after solvation of a dry geometry. In 
our computations, we estimated the water volume in a few realizations of each 
configuration of interest, several measures of p n were collected and used to 
compute an average value and the corresponding s.d. The latter s.d. generates 
uncertainties when computing the volume V w , and consequently uncertainties of 
the scaling parameter (horizontal error bars in Fig. 3). Alternative approaches 
based, for instance, on Monte-Carlo integration can also be used to calculate V w 
accurately. 
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Characterization of nanoconstructs and NPs. Details on the fabrication and 
characterization of the discoidal SiMPs and the SPIOs are provided in the 
Supplementary Methods and Supplementary Figs 33 and 34, together with the 
protocols for the loading of SPIOs into SiMPs. 

Relaxivity measurements. In vitro relaxation times were measured in a Bruker 
Minispec (mq 60) bench-top relaxometer operating at 60 MHz and 37 °C. The 
transverse (T 2 ) relaxation times were measured using the Carr-Purcell-Meiboom- 
Gill sequence. 

Analysis of relaxivity enhancement. The function / EDX (X,Y) of the two spatial 
coordinates X and Y is built from the EDX image for iron (Fig. 5c) and denotes the 
position of the centre of a generic image pixel. Let Nspio be the average number of 
SPIOs within a single SiMP. Hence, a surface density function at a pixel scale for 
the SPIOs can be written as 



Pareal- 



/edx(^, Y) Nspio 
J / EDX (X, Y)dXdY 



(12) 



where A pix and A SiMP are the surface area of a single pixel and the surface area of a 
SiMP, respectively. The function p area i has been reported in the (left) bottom part of 
Supplementary Fig. 35. SiMP pores that are loaded by SPIOs can be automatically 
recognized by finding the local density map of the EDX non-zero elements (p p0 int)- 
The density function in equation (12) depends on the image pixel size (that is, a 
quantity that is not representative of the phenomenon under study), whereas a 
more meaningful function should be based on the pore characteristic dimension 
instead (representative of the confinement length). Hence, the following surface 
density function at a pore scale for the SPIOs is introduced: /? SPIO (X,Y) = 
Cp V oint(X,Y), with C being a constant to be determined. To this end, owing to mass 
conservation, it is imposed 



J p SPI0 (X, Y)dXdY= J p areal (X, Y)dXdY, 



hence, 



J> areal (X, Y)dXdY 
7p point (X, Y)dXdY' 



(13) 



(14) 



with A p being the surface area of a representative SiMP pore (see also the (right) 
bottom part of Supplementary Fig. 35). The total number of SPIOs within the 
representative pore can be computed as N pore = f A p SPIO (X, Y)dXdY. Hence, the 
water volume near the SPIOs (for the latter pore) V w can be estimated as 
V w =A p hem — Vspio-Npore> where h is the SiMP height, while 0< ca\< 1 indicates 
the fraction of the SiMP pore that has been effectively occupied by SPIOs. Clearly, 
for homogeneously distributed SPIOs within the pore, c m = 1. For the considered 
representative pore, an estimate of the average scaling parameter 9 can be obtained 
as: 9=1 -exp( -0')with 



-, Np 0r e S SP io ^SPIO + 71 Dh Cffli <S Si 



(15) 



where D is the average diameter of the pores. Finally, a map of the scaling para- 
meter 9 = 9(X,Y) is derived by the map Pspio as G(X y Y) = C Pspio(X,Y), where the 
constant C is determined by imposing 



— / 0(X, Y)dXdY=9, 



hence 



C'= 



9A p 

/Pspio(*> Y)dXdY' 



(16) 



(17) 



On the estimate of the scaling parameter map 9 = 9{X,Y), the computation of a 
corresponding map for the diffusion coefficient D = D(X,Y) is straightforwardly 
achieved by adopting the law suggested in this work. 

The above procedure is slightly sensitive to the choice of the representative 
pore. However, from our computations no significant changes in the final result 
(map of 9 = 9(X,Y)) was experienced by making difference choices. In 
Supplementary Fig. 36, we report both the map of the scaling parameter 9 = 9{X,Y) 
and the corresponding D = D(X,Y) at different pore filling. 

Let us consider two contrast agents with the same iron molarity M Fe . Let the 
first one be based on SPIOs homogenously dispersed in bulk water, while the 
second one adopts the same SPIOs (nanoconfined) within SiMPs. Given the 
definition of r 2 (Supplementary equation (9)), we can define a relaxivity 



enhancement as follows: 



En~ 



Since 1/T 2 > 1/T 0 - 



(r 2 ) 2 1 






M Fe 


(rO, 


M Fe | 




',-ifc 



(18) 



En 



(19) 



Finally, using the outer sphere theory (see Supplementary Discussion), the 
enhancement can be recast as 



En 



UJ 2 = 9 tJ 2 fc(Aco r ) = v 2 Di 

(^) fa i-S^f S»r 

Consistently with the above procedure, we aim at estimating the map En(X,Y): 

gi » 2 (X, Y) 
oi D 2 (X, 7)' 



En 



(20) 



(21) 



where D x and V\ are the self- diffusion coefficient of bulk water and the average 
volume fraction of SPIOs, respectively. Let us assume that u 2 (X, Y)=C" p mo (X, Y) 
with j- J u 2 (X, Y)dXdY=u; hence, 



c // = » a p 

/Pspioft Y)dXdY' 



(22) 



and 



(23) 



_ Npore^SPIO 

A h Cfiu 

In Fig. 5d we report the function En(X,Y) corresponding to the SiMP particle of 
Fig. 5c. 

Finally, the average enhancement due to the entire SiMP can be estimated as 



En= 



f smo {X, Y)= 



(24) 



(25) 



J En(X,Y)f smo (X,Y)dXdY, 

AsiMP 

where f S pio(X,Y) is the following distribution function 

PsPiofe Y) 
' I P mo (X,Y)dXdY' 

For the case in Fig. 5d, we find En ~ 2.7. In Fig. 5c, we consider SiMPs filled by 
5 nm SPIOs, the transverse relaxivity of which was experimentally assessed as 
(r 2 ) 2 = 270 ± 73 mM _ 1 s _ l . Moreover, 5 nm SPIOs ( Sigma- Aldrich) show trans- 
verse relaxivity in bulk conditions of water (r 2 )i = 107 ± 24 mM ~ 1 s ~ 1 . As a result, 
the experimentally measured enhancement for the agent in Fig. 5c (compared with 
SPIOs in bulk water) is En w 2.52, which is in excellent agreement with the above 
prediction, fully based on the EDX signal (iron) and our model. 
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